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Abstract 

For a class of coupled limit cycle oscillators, we give a condition on a linear coupling operator 
VO that is necessary and sufficient for exponential stability of the synchronous solution. We show that 

with certain modifications our method of analysis applies to networks with partial, time-dependent, and 

, , nonlinear coupling schemes, as well as to ensembles of local systems with nonperiodic attractors. We 

also study robustness of synchrony to noise. To this end, we analytically estimate the degree of coherence 
of the network oscillations in the presence of noise. Our estimate of coherence highlights the main 
ingredients of stochastic stability of the synchronous regime. In particular, it quantifies the contribution 
of the network topology. The estimate of coherence for the randomly perturbed network can be used 
as means for analytic inference of degree of stability of the synchronous solution of the unperturbed 
i i deterministic network. Furthermore, we show that in large networks, the effects of noise on the dynamics 

of each oscillator can be effectively controlled by varying the strength of coupling, which provides a 
powerful mechanism of denoising. This suggests that the organization of oscillators in a coupled network 
may play an important role in maintaining robust oscillations in random environment. The analysis is 
complemented with the results of numerical simulations of a neuronal network. 
PACS: 05.45.Xt, 05.40.Ca 
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x t = f(x t ) + cP(t)w t , x : R 1 R n , (1.1) 

l _> 



On 1 Introduction 

O 

Consider a dynamical system forced by small noise 

>< 
U 

where function f : W 1 — > W 1 is continuous together with partial derivatives up to the second order, P : 
W 1 is a bounded continuous function of time, w t is a n— dimensional white noise process, and small a > 
is the noise intensity. We call ( |1.1[ ) a local system and denote ( l.l[ )o the underlying deterministic system 
( |1.1[ ) with a = 0. For a > 0, stochastic differential equation (JTTT ) is understood in the Ito's sense [1J. Many 
models of (bio)physical phenomena are formulated as the coupled networks of N identical local systems 
<0>: 

x t = f(x t ) + gDx t + <rP(t)w u (1.2) 
where x = (xW,x^\ . . . ,xW) €R Nn ,f(x) = (f (x^)), f (x( 2 )), . . . , f(x( N ))) G M Nn and 
P{t) = In <8> P(t), and wt = ( , , . . . , ) , are independent copies of n— dimensional Brow- 
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nian motion. The coupling is implemented by a linear operator D : 'K Nn — > ~R Nn and g > is inter- 
preted as the strength of coupling. D may depend on time. Our only assumption on D is that it leaves 
the diagonal invariant, i.e., if x = £(t) solves ( | 1 - 1 1 ) q then x = £(t) := In <8> £(t) solves (1.2)o. Here, 
1 N = (1,1,..., 1) T g 1^ and (g> denotes the Kronecker product, so that £(t) = (£(t), £(t), . . . , £(t)) is 
a solution of the coupled system. Such solutions, when asymptotically stable, feature synchronization, an 
important mode of collective behavior. Suppose x = £(t) is an asymptotically stable solution of the local 
system ( |1 . 1 [ )q. Under what conditions on the coupling operator D, x = In ® £(t) is an asymptotically 
stable solution of ( 1.2 )o? What information about D is important for synchronization properties of the cou- 
pled system? What determines the rate of attraction of the coupled limit cycle and its robustness to noise? 
Clearly, the answers to these questions depend in a nontrivial fashion on the properties of the local systems, 
network topology, and the type of interactions between local systems. 

Many different approaches have been proposed to study synchronization: asymptotic analysis EJ [3], 
phase reduction 0|5]], constructing Lyapunov functions iflOl and estimating Lyapunov exponents 0, invari- 
ant manifold theory SHI, and graph-theoretic techniques [9]. Above we selected just a few representative 
studies illustrating these approaches. For more background and more complete bibliography, we refer an 
interested reader to recent monographs tiTTl[T2l[T6l . For weakly coupled networks, i.e., when < j < 1, 
there are effective perturbation techniques for studying synchronization such as asymptotic approximation 
of the Poincare map and the method of averaging. These ideas underlie widely used method of phase re- 
sponse curves and Kuramoto's phase reduction Bill. For a review of techniques available for studying 
synchronization in weakly coupled networks we refer to Chapter 10 in |[T3l and references therein. When 
the coupling is moderate g = 0(1) or strong, these methods do not apply. Moreover, the mechanisms 
for synchronization in weakly and strongly coupled networks are different. In the case of strong coupling, 
a prevalent approach for studying synchronization is to use the properties of specific (albeit important in 
applications) coupling schemes such as global all-to-all coupling (see, e.g., |[T4l "). local nearest-neighbor 
coupling, or more generally schemes resulting from discretization of Laplace operator |2l|6]. For these 
network topologies, one can use explicit information about the spectra of the coupling matrices; in addition, 
the former scheme has strong symmetry properties that can be used in understanding network dynamics. 
Networks with general coupling operators have been studied by Pecora and Caroll [7 ] and by V. Belykh, 
I. Belykh, and Hasler SHU. The master stability function, constructed in [7J, uses spectral properties of a 
given coupling operator to determine whether the synchronous solution is stable. Practical implementation 
of this method relies on numerical computation of matrix eigenvalues. Analytical sufficient conditions for 
synchronization derived in |£l [151 use graph theoretic interpretation of the coupling operator to construct 
Lyapunov functions controlling the growth of perturbations of the synchronous solution. In this Letter, we 
look for an analytic or rather algebraic description of coupling operators that endow synchronous solutions 
with exponential stability. We define a class of matrices, dissipative matrices (cf. Definition [TJ, and show 
that these matrices generate exponentially stable synchronous solutions once the coupling strength exceeds 
a certain value. A random dissipative matrix shown in Section [4] suggests that many dissipative matrices 
do not fall into the class of coupling operators analyzed in (HQS. Therefore, by identifying dissipative 
matrices we have substantially extended existing knowledge of linear coupling operators that enforce syn- 
chrony in coupled networks. Surprisingly, dissipative matrices admit an explicit algebraic characterization: 
Theorem [2] relates all dissipative matrices to a discrete Laplacian, justifying common interpretation of the 
Laplacian as a prototypical diffusive coupling. To highlight the main ingredients of our treatment of syn- 
chronization, in Section [2] we present the analysis in the simplest meaningful setting when the coupling is 
linear and stationary. In Section 3, we discuss how to apply our method to problems with time-dependent 
and nonlinear coupling operators, and ,importantly, to networks composed of local systems with nonperiodic 
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attractors. 

Adequate description of many physical phenomena requires including stochastic terms into differential 
equation models. In the context of synchronization this leads to an important question of robustness of 
synchrony to noise. This is the second problem investigated in this Letter. We analytically estimate the 
coherence of the coupled oscillations in the presence of noise. The estimate is tight. It reflects the main 
ingredients of robustness of synchronous oscillations to noise. In particular, it quantifies the contribution 
of the network topology to the stability of the synchronous solution. As a related result, we show that in 
large networks the effects of noise on oscillations can be reduced substantially by increasing the strength of 
coupling. For networks of simpler elements, so-called integrate and fire neurons, the denoising property was 
shown in [19]. The present Letter extends the result of fl9l to systems of coupled limit cycle oscillators. 
Finally, in Section 4 we illustrate our findings with a discussion of the dynamics of a concrete biophysical 
model, an ensemble of neural oscillators. 

2 The analysis 

We start by specifying the structure of the coupling operator. We call coupling operator D separable if 

d D l. i) < :> Vx v .l c :s" x ". an 

Matrices D and L play distinct roles in the network organization. D reflects the global architecture: what 
local system is connected to what. L specifies how the coupling is organized on the level of a local system: 
roughly, what local variables are engaged in coupling. The separable structure of the coupling is important. 
It translates naturally to the stability analysis of the synchronous state. The condition that the diagonal is 
invariant for separable coupling translates to In G ker D. Moreover, if the network is connected then 
ker D is one-dimensional. Thus, we are led to the following condition 

D G K = {M G R NxN : ker M = Span {1 N }} • (2.2) 

Further, we assume that L is symmetric semipositive definite, i.e., L T = L and x T Lx > Vx G M n . D is not 
assumed symmetric. For simplicity, we keep D and L constant. In Section [3] we explain our results for 
time-dependent and nonlinear coupling. 

If L is positive definite, we say that the coupling is full (rank), otherwise we call it partial (rank). The 



distinction between the full and partial coupling is important for synchronization properties of (1.2). The 
following examples illustrate how full and partial coupling arise in physical models. Let 



/ -JV + 1 1 1 ... 1 \ 
1 -JV + 1 1 ... 1 

V i i i ... -jv + i J 



(2.3) 



I be the n x n identity matrix, and I' = diag(l, 0, 0, . . . , 0) G M nxn . D := Dj ® I in ( 1.2) implements 
all-to-all coupling through all variables, while D := Di (g> I' engages only the first variables of the local 
systems. The former is an example of the full coupling, whereas the latter is that of the partial coupling. The 
stability analysis of partially coupled systems has to deal with the degeneracy of the coupling matrix D (due 



to the zero eigenvalues of L). A reader wishing to gain better physical intuition for coupled system ( 1.2 1 and 
( |2.1| ) before embarking on analysis, is referred to the discussion of a compartmental model of a neuron in 
Section H 
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To study synchronization in ( 1.2l-(2.2), we derive the equation for the phase variables of the coupled 
system. To set up the notation, we review phase reduction for a single oscillator. Let x = £(t) denote 
a periodic solution of ( |l.l| )o with the least period 1. By O = {x = £(t), t G S 1 = M/Z} we denote the 
corresponding orbit. Along O we introduce an orthonormal moving coordinate frame (cf. iTTTl ): 



(2.4) 



where the first vector is a unit vector moving along O, v(9) = £(0) £(9) 



-l 



The change of variables 



x = m + Z(0)p, 2(9) = col(zi(fl), . . . , z n _i(0)). 



(2.5) 



defines a smooth transformation x h-» (6*, p) G S 1 x M n_1 in a sufficiently small neighborhood of O (cf. 
Theorem VI. 1.1 in IPT71 ). Using Ito's formula, in the vicinity of the periodic orbit, we rewrite ( |1.1[ ) in 
terms of 9 and p (cf. (2.5)) and project the resultant equation onto the subspaces spanned by v(0) and 
{zi(0),...,z n _i(0)} to obtain 



9\ = l + (7hl(0 t )P(t)w t + ..., 

Pt = A(0 t )p t + ah 2 (0 t )P(t)w t + 



(2.6) 
(2.7) 



where 



A(9) = Z(9) T 



-82(9) 
89 



+ Df(£(0))Z(0) 



(2.8) 



h 1 (0)=v T (0)|e(0))|- 1 , h 2 (0) = Z T (0). 



The detailed derivation of the system of equations near a periodic orbit of a deterministic system in the 
moving coordinates can be found in the proof of Theorem VI. 1.2 in IfTTl . The treatment of the stochastic 
case requires Ito's formula, which does not affect the leading order terms written out in (2.6 1 and (2.7 1 (see 
the proof of Lemma 4.1 in |fT8l for details). The solution of an initial value problem for (2.6 1 and (2.7i 
yields a real- valued function 9 t . The phase of the oscillations is given by (9 t mod 1) G S 1 . It will be more 
convenient to work with 9 t , which we will call a phase variable, rather than with it's projection on S 1 . 
Assume that the eigenvalues of A s (9) = A(9) + A T (9), \i(9), i = 1, 2, . . . , n — 1, are negative 



maxAi(6») < -A < 0. 
ees 1 



(2.9) 



By applying the phase reduction to each oscillator in the network, in complete analogy to (2.6), we derive 
the phase equations for the coupled system 

i(0 



C =l + (7hl(0 t W )P(t)w{ 



(2.10) 



where dij denote the entries of D (cf. (2.1 )). The expression for the coupling terms in (2.10) simplifies to 



v T (0( l )) | 



£(0G))_ e( 0(i)) 



g(g( | )) T U(g (i 

i(9&) 2 



4 



(flCi) _ 0(0) + . . . 



(2.11) 



Here and below, we ignore quadratic terms O ((6>W — 9^) 2 ). By plugging (2.11 1 in (2.10l, we arrive at the 
following system of equations 



B t = In + gV N (9 t )T>9 t + aH 1 {9 t )P{t)w t + ..., 



(2.12) 



where 



and 



V N (9) = diag (l(0W), l(6&), . . . , 1(9^)) , H x {9) = diag (h^ 1 )), . . . , 



(N)- 
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1(9) := i(9) m T me)=v(9) T Lv(9). 
Next, we derive the system for the vector of the phase differences 



S9 



(N-l) 



where (N — 1) x N matrix S is defined by 



/-ll 
-1 1 



(i) _ fl (i+l) _ fl (i) 



\ 




(2.13) 



(2.14) 



(2.15) 



V ... -1 1 / 
Multiply both sides of ( |2.12| ) by S and note that 

S1 N = and SV N # = V N _iS# + 0(\<f>\), 



where is defined in (2.14). From (2.12) and (2.16) we have 

4>t = gV N ^(9 t )t)<j)t + crSH l {9t)P{t)wt + . 



(2.16) 



(2.17) 



wh ere D is defined by SD = DS. For D G /C, D is well-defi ned (c f. Appendix [19]). In the derivation 
of (2.17), we treat ed ter ms ~ 1 0| = O(|0| 2 ) inherited from (2.16) as higher order terms. Since #W = 
gW+Q(M) (cf. gig), 

Vjv-i(e) = /(^)I N -i + O(|0|), 
#i(0) = In + O(|0|), 



and ( 2. 17 1 is reduced to 



<?/(^ (1) )D0 t + <tS(I n ® /ii(^ (1) ))(I N 8) P(t))w t + ... 
gl^Wt + (7S(I N ® /ii(^ (1) )P(t))«; t + . . . . 



(2.18) 



In (2.18), O(|0| 2 ) terms are treated as higher order, because they do not affect exponential stability of the 
synchronous solution. Note how separable coupling translates to the structure of the phase equation. Matrix 
D, which carries the information about the network topology effectively determines the stability of the 
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synchronous solution. Semipositive matrix L enters the factor 1(9^) (cf. (2.13l). The stability of (2.18) is 
determined from the homogeneous deterministic system: 



j H = gi(el 1) )i) 



(2.19) 



where by 0^ we mean the first component of the solution of deterministic equation < |1.2[ )o. We continue our 
analysis assuming that L > 0, i.e., the coupling is full. In Section [3] we comment on how our results apply 
to partially coupled systems. Thus, l(o[^) > a > and after changing the independent variable we have 

4> T = gb<t> T . (2.20) 

For exponential stability of synchronous solution, the symmetric part of D must be negative definite. This 
motivates the following definition. 



Definition 1. Matrices from 



T> = |m G K. : x T Mx < Vx G R jV-1 /{0}} 



(2.21) 



are called dissipative. 



Thus, we arrive at the first conclusion of this Letter: synchronous solution of ( 1.2)o-(2.2i is exponentially 
stable iff D is dissipative. When studying (1.2)q-(2.2| it is tempting to relate the stability of synchronous 
solution to the spectrum of D. It is important to realize that it is the spectrum of D that is responsible for 
synchronization. Remarkably, dissipative matrices admit an explicit characterization. 



Theorem 2. D G V iff 

D = QA , A = S T S 

for some Q G M. NxN with negative definite symmetric part. 



(2.22) 



Proof. Suppose {2.22) holds. Let D = SQS T . Then 

SD = SQS T S = DS. 

Furthermore, for any x G R Ar_1 /{0} 

(Dx,x) = (SQS T x,x) = (Qy,y), y = S T x, 



(2.23) 



where by (-, •) we denote the inner product in E^ -1 . We will use the same notation for the inner product in 
any other Euclidean space used in this Letter, e.g., R N and C^. Note that y = S T x ^ 0, because the rows 
of S are linearly independent. Thus, from ( |2.23| ) we conclude 

(Dx,x) < Vx G R^VlO}, 

i.e., D G V. 

Conversely, suppose D G V. Note that on the orthogonal complement of ker Ao, (ker Ao) _ 
Aq is invertible and define Q G M. NxN as follows: for x G ~M> N let 



Qx 



D(A | lN x)- 1 x, xGIn 1 , 

—x, x G Span{l]\j}- 



In 1 , 



(2.24) 
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We show that the symmetric part of Q | 1n ± (and, therefore, Q itself) is negative definite. For any z £ 
In^/IO} there exists x £ M 7V_1 /{0} such that z = Aox, because Ao \ 1 ± is invertible. Moreover, such 
x, can be chosen from 1n j -/{0} because ker Ao = Span {In}- Thus, 

(Qz, z) = (D(A | 1n x ) _1 A x, A x) = (Dx, S T Sx) = (SDx, Sx) = (DSx, Sx) < 0. (2.25) 

Here, we used the fact that R(Ao) = In" 1 and Sx 7^ (because x £ (ker S) -1 ). The combination of (2.24 1 
and ( |2.25| > yields ( |2.22| ). 

□ 

Theorem [2] gives explicit and for separable coupling exhaustive characterization of coupling matrices 
that generate exponentially stable synchronous solutions. Synchronization is often attributed to systems 
with diffusive coupling that are obtained by discretizing elliptic differential operators or, more generally, 
differential operators modeling diffusion on graphs. In this respect, it is remarkable that Theorem [2] relates 
all dissipative matrices to the discrete Laplacian 



/ 1 



An 



S T S 



■1 



V 



2 









\ 





(2.26) 



-1 1 J 



This is consistent with the common interpretation of the Laplacian as a prototypical elliptic operator. The 
explicit characterization of dissipative matrices by (2.22 1 is the second result of this Letter. Before turning 
to the question of the robustness to noise we state two corollaries of Theorem [2] The first corollary gives 
a convenient computational formula for D, while second one characterizes the spectrum of a dissipative 
matrix. 



Corollary 3. For D £ £>, D = SQS T , where Q satisfies \2.22\ . 

Corollary 4. If D € T> then D has a zero eigenvalue of multiplicity 1. All nonzero eigenvalues of D are 
real and negative. 



Proof. Matrix D has a simple zero eigenvalue by (2.2). Suppose A £ C is a nonzero eigenvalue of D and 
x £ C n be a corresponding eigenvector 

QA x = Ax. (2.27) 



Note that x ^ Span{l_N}. By multiplying both sides of (2.27 1 by Aox 7^ 0, we have 

(QA x, A x) = A(x, A x) = A(Sx, Sx). 



(2.28) 



The scalar product multiplying A on the right hand side of (2.28 1 is positive, while (QAox, Aox) is negative. 
Therefore, A < 0. □ 

Having understood the mechanism for synchronization in the deterministic network (1.2)n. We now 
turn to the question of robustness of the synchronous regime to noise. To this end, we return to (2.18 1. For 
the remainder of this Letter, D £ V. For small a > 0, on a finite time interval solution of (218 1 can be 
expanded as 

4>t = 4>t + afa + ■ ■ ■ , (229) 



where deter minis tic function <pt solves (2.19) (cf. Theorem 2.2, [20]). Since <fi = is an exponentially stable 
solution of (2.19), we take 4> t = 0. The leading order correction a(p t is a Gaussian process, which for small 
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a > approximates <pt on a finite interval of time. This specifies the scope of applicability of our analysis. 
In particular, we are not concerned with large deviation type effects which become relevant on much longer 
timescales. 



From ( 2. 1 2 1 we have 



9^=t + 0(a). 



By plugging ( [2T29] ) and p30| ) in ( pU8j ), we have 

4> t = gl(t)t)4> t + S(I N ® ^i(t)P(t))?ii + 



After changing time to r = J* Q * l(s)ds and ignoring higher order terms, (2.32 1 is rewritten as 
$ T = gt)4> T + S(I N <g> /n(r)P(r))«; T , fc(r) := P(r) := P(t(r)). 



(2.30) 



(2.31) 



(2.32) 



The solution of (2.32 1 subject to initial condition 4>q = is a Gaussian random process with zero mean and 
covariance matrix (cf. §5.6, HI) 

cov o T I e 9(T " s)6 S(I N ®Ms)P(s))(I N ® h{s)P(s)) T S T e^ T - s)f)T ds 

h 2 (s)e 9{T - s)ti Ae^ T - s)f)T ds, (2.33) 

where A := SS T and h 2 {t) = In <8> (h(t)P(t)P(t) T h T (t)) is a nonnegative scalar function. Using standard 
properties of trace and the mean value theorem, from ( |2.33 1 we have 

TV COV 6-r = 



h 2 {u)Tr {Ae^'^jdu 
/i 2 (C(T))Tr{Ae 9rf)S e~ 9ui)a du} 



h 



! (C(r))Tr {-Ag-^f)*)- 1 



In-i - e 



9 rD s 



M (r)^ + 0(e- c -) 
5 



(2.34) 



where continuous function £(r) is due to the application of the mean value theorem, /x(t) := /r(£(r)), 
D s := D + D T , and 

(2.35) 



«(D) = Tr{-A(D s )- 1 } 



Note that < /i(r) < M is uniformly bounded. Define the average variance of the variables <p[ k \ k = 
1,2 JVas 

1 7V_1 1 

Yl var ^ = A 73I Tr cov </>*■ (236) 
fe=i 



var <pt 



N 



Equation (2.34) yields an important estimate for the network variability 



var (f> T m cr 2 var (f> T = <j 2 //(t) 



«(D) 



(2.37) 
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where c\ is a positive constant. Nonnegative function /x(r) reflects the properties of the local system such 
as geometric properties of the limit cycle and matrix P multiplying stochastic term, whereas k(D) captures 
network topology. For a network of fixed size, vaf 4>t can be made arbitrarily small by taking large g. 
Moreover, by Chebyshev's inequality, for any S > 0, 



ct 2 MNk(D) 



as g 



oo, 



i.e., for strong coupling the phases of individual oscillators can be localized within arbitrarily narrow bounds. 
The control of the coherence by varying the coupling strength is more effective in networks with smaller 



k(D). Thus, (2.37 1 shows explicitly the factors controlling the coherence in the presence of noise. Moreover, 
«(D) quantifies the contribution of the network topology to the stability of the synchronous state. This is 
the third conclusion of this Letter. 

What features of the network topology are captured by k(D)? We first go over the ingredients of the 



formula for k(D) (cf. (2.35 )). Matrix A is a Laplacian: 



SS J 



2 


-1 


. . 


. 





-1 


2 


-1 . . 


. 











. . 


. -1 


2 



(2.38) 



Unlike Ao, A is nonsingular. The following examples show that «(D) can change by many times for 
networks with different topologies. Let Di be as in (2.3) and D2 = — Ao (cf. ( 2.26| >). Di and D2 are 
coupling matrices corresponding to the graphs modeling all-to-all and nearest-neighbor interactions in the 
network. By direct verification, 



-2N I N _i and D| = -2A. 



(2.39) 



By plugging the explicit expressions (2.39 1 in (2.35 1, we find 



«(Dj.) = 1 + 0(N~ L ) and k(D 2 



N- 1 



(2.40) 



Note that the all-to-all topology features a significant reduction in k compared to the nearest-neighbor cou- 
pling. This reduction is proportional to the ratio of the degrees of the corresponding graphs: 2 - for the 
nearest-neighbor and (N — 1) - for the all-to-all coupling. Thus, (2.37 1 suggests that networks with the 
higher density of connections are more robust to noise. 



Equation (2.37 ) estimates the variability of the phase differences, revealing the main factors contributing 



to robustness of synchrony to noise. If D is symmetric, Equation (2.37) can also be used to estimate the 



variability of the phase variables 0^ themselves. We follow the method used in |T9l for a related problem. 
First, we derive the equation for the average phase of the coupled system 

e = n T 9, v = N- 1 i^. 



By multiplying both sides of (2.121 by 77 , we have 



1 N 



t = l + <7h 1 (t)P(t)X t + ..., 



(2.41) 



i=i 
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Here, we used the following approximations 

V N {9) = Z(0 t (1) )I N + O(|0|) and t (1) = t + 0(a). 

As a linear combination of independent Gaussian processes , X(t) is distributed as iV _1 / 2 Wt, where w t 
is a n— dimensional Brownian motion. Thus, 

|t = 1 + ah l( t)P(t) Wt+ (2A2) 

and 



N 



var 9 t ^ jf h 1 (s)P(s)P(s) T h 1 (s) T d S < t G [0,T], (2.43) 

where C\ > does not depend on iV and T. Next, by noting as in [19] that each phase variable #W can be 
represented as a linear comb inatio n of t he average phase 9 and phase differences <p^\ 4>^ 2 \ ■ ■ ■ , 4>^ N \ we 



estimate var #W using (2.43 1 and (2.37 1. Omitting further details, which are the same as in step 3. of the 



proof of Theorem 3.1 in 11191 . we state the final result 

„m •> (C 2 T C 3 Nk(D)\ 
max sup var t (lj < a 2 + — ^ , (2.44) 

» te[o,T] ' V N 9 J 

where C 2 and C3 are positive constants independent from N, g and T. The first term on the right hand side 



of (2.43 1 can be made arbitrarily small by increasing N, while the second term decreases for increasing g. 



Therefore, in large networks, the effects of noise on oscillations can be controlled by varying the strength of 



coupling. The two terms on the right hand side of (2.43 ) represent two main ingredients of the mechanism 



of denoising: the first term is due to the averaging of statistically independent stochastic forces acting on 
connected local systems, whereas the second term reflects the dissipativity of the coupling. The latter is 
a critical property of the coupling operator that underlies both synchronization and denoising in coupled 
networks. 

The analysis of the phase equations above produced a necessary and sufficient condition for synchro- 
nization in systems with separable coupling and gave a compact explicit estimate for the spread of phases 
of coupled stochastic oscillators. For the phase equations to be valid, the trajectory of the coupled system 
must remain close to the limit cycle. To complete the analysis, we consider the system of equations for 
Pt = (fit 5 Pt 1 ■ ■ ■ 1 Pt^)- The derivation of the system for p is completely analogous to that for B t (cf. 
( |2.12| l). We omit the details and state the final result 



Pt = [A{t) + gD(t)] Pt + ah 2 (t)P(t)w t + ..., (2.45) 



where A(t) = I N ® A(t), D(t) = D <g> (Z T (t)LZ(t)). Matrices Z(t) and A(t) are defined in (2.5 1 and 



(|2.8[> respectively. By Vazhevski's inequality Ell , the combination of (|2.9[>, D G V, and L > implies 



exponential stability of the equilibrium at p = in (2.45 )q. Therefore, on finite time intervals with over- 



whelming probability, \p t \ remains small, provided |po| and a > are sufficiently small. This justifies the 
phase reduction that we used above. Note that this conclusion holds for both full and partial coupling. 

3 Generalizations 

The analysis of this Letter admits several important generalizations. 
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A) Partial coupling. If the coupling is partial, nonnegative function /(•) in (2.32 1 in general takes zero 
values. Dealing with the degeneracies in ( 2.32| ) requires additional care. For a common in applications case 
when /(•) has isolated zeros, with technical modifications one can get a qualitatively similar estimate to 
((237). 

B) Time-dependent coupling. Our analysis remains unchanged if instead constant L one uses a bounded 
measurable function of time. The definition of the full coupling is then modified to x T L(t)x > ax T x for 
some a > and Vx G M n /{0} uniformly in t > 0. Likewise, D can be taken time-dependent as long as 
D(t) G T> for all t. In this case, exponential stability of <p = follows from (2.20 1 if we require that all 
eigenvalues of D s (i) = D(i) + D T (i) are negative and bounded away from zero uniformly in t: 

x T D s (t)x = x T D(t)x < - 7 x T x Vx G R JV_1 /{0}, t > 

for some 7 > 0. By Theorem]!} such matrices can be written as D = Q(t)Ao, where Q(t) is such that 

x T Q(t)x < -7x T x Vx G M^/jO}, * > 

for some 7 > 0. Also, in the time-dependent case, one can get a slightly weaker but qualitatively similar 
estimate on var 4>t (cf- (237 1). 

C) Nonlinear coupling. The analysis can be extended to the systems with nonlinear coupling of the 
following form 

N 

x« = f(x«) + g £ d u (t)L(x«), i = 1, 2, . . . , N, (3.1) 



where D(t) = (dy(t)) G K for every t > and L : W l -)• R n is a smooth function. Since D(t) G K, ( 3.1) ) 
can be rewritten as 

N 

x« = f(x«) + g dij(t) (C(x«) - L(x«)) , i = 1, 2, . . . , N. (3.2) 

0=1 



Suppose x = ^(t) is a periodic solution of the local system x = f(x). Then x p = In ® solves (3.2 1. In 
a neighborhood of x p , (3.2) can be rewritten using Taylor's formula 



N 

x« = f(x«) + g dy(t)L(t)(x® - xW) + 0(max |x^ - x^] 2 ), i = 1, 2, . . . , N, 
3=1 kJ 



or, equivalently, 

± = f(x) + <?(D(t) (8) L(t))x + . . . , where L(t) := (3.3) 

Equation ( |3.3| ) is now in the form, for which the analysis of Section 2 applies (cf. B) above). 

D) Nonperiodic attractors. Moving coordinate systems similar to the one used in this Letter can be 
introduced in the vicinity of locally invariant sets of more general nature. For example, motions along 
certain attracting normally hyperbolic slow manifolds admit a similar description (cf. |[T7ll23l ). Thus, our 
analysis can be adopted to study synchronization in a more general setting. Furthermore, in the case when 
the attractor of the local system is periodic and synchronization takes place, the analysis of Section 2 yields 
a precise description of the asymptotic behavior of trajectories of the coupled system. Specifically, the 
attractor of the coupled system includes a periodic orbit that is a direct product of the periodic orbits of the 
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local systems. If one is only concerned with synchronization and does not aim at describing the asymptotic 
behavior of the coupled system, then the transformation to moving coordinates is not needed. To outline the 
analysis for this case, let x = £(t) be a solution of the local system. This solution is not necessarily periodic. 
Instead, we assume that x = £(t) does not leave a bounded domain in R n . Then £(t) = In <8> £(t) solves 
the coupled system 

x = f{x)+g(-D^L)x,x = (x^\x^,...,x^)eR Nn . (3.4) 
Linearization about x = yields 

i = (I N ®A(t) + p(D®L))or + ..., A(t) = ^|^. (3.5) 

ax 



By multiplying both sides of <\3.5) by S = S (g) I, we get the equation ofy = Sx: 

y = (I N _i<g)A(t)+s(D<g)L))2/ + .... (3.6) 



Asymptotic stability of y = implies synchronization for (3.4i. A sufficient condition for asymptotic 



stability of the trivial solution of (3.6 1 is that the eigenvalues of symmetric matrix 

B = I N _i ® A s (t) + gf) s <g) L (3.7) 

are negative and bounded from zero uniformly in t > 0. Since |A s (t)| is bounded, the desired property for B 
for large g follows from D 6 V and L being symmetric positive definite. Thus, we get a sufficient condition 
for synchronization for the full coupling. In the partial coupling case, D s ® L has N x (n — rank (L)) zero 
eigenvalues and one has to make sure that they do not give rise to negative eigenvalues of B. The condition 
for this can be obtained from the well-known formulas for the perturbations of the eigenvalues of symmetric 
matrices (cf. Appendix in [24]). A more complete analysis of synchronization for the partial coupling case 
will be given elsewhere ESI . 

4 Numerical example 

To illustrate the analytical results of this Letter, we use a nondimensional model of a pacemaker neuron from 
(221: 

eiJ® = g 1 (vP)(E 1 -vf ) )+g 2 (u^)[E 2 -v^)+g 3 [E 3 -v^)+I^+a4 i \ (4.1) 
u? = "Liv^fa-v^-^J+I®, ^ = 1,2,...,N. (4.2) 

Dynamic variables «W and ttW represent membrane potential and calcium concentration in a given compart- 
ment of an axon or a dendrite of a neural cell. The compartments are sufficiently small so that the membrane 
potential and calcium concentration can be assumed constant throughout one compartment (Fig.[T^). Terms 



on the right hand side of the voltage equation (4.1 1 model ionic currents: a calcium current, a calcium de- 



pendent potassium current, and a small leak current. In addition, small white noise is added to account for 



random synaptic input or other fluctuations. The equation for calcium concentration (4.2 1 takes into account 



calcium current and calcium efflux due to calcium pump. The ionic conductances are sigmoid functions of 
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cm i □ ) 

if 



i) u d) W ) 



v (N)(N) 




Figure 1: Schematic representation of spatial structure of a compartmental model: (a) linear cable, (b) 
branched cable. Dynamic variables and r', i = 1, 2, . . . , JV, approximate voltage and calcium concen- 
tration in each compartment. 

the voltage and calcium concentration 

gi {v) = |(l + tanh(^)), (4.3) 

92{u) = (4.4) 

We briefly comment on the meaning of the model parameters: #1,2,3 and £4,2,3 stand for maximal conduc- 
tances and reversal potentials of the corresponding ionic currents; 01,2,3 are constants used in the descrip- 
tions of activation of calcium and calcium dependent currents; ui and e are certain constants that come up 
in the process of nondimesionalization of the conductance based model of a dopamine neuron (see [22J for 
details). The values of the parameters that were used in our simulations are given in the appendix to this 
Letter. 

The coupling terms 

JV 

J« = g^d^-vP), (4.5) 

i=o 

J« = 6j2diM j) -4 ] ) (4-6) 
3=0 

model electrical current and calcium diffusion between adjacent compartments respectively. In case of 
a linear cable geometry of an axon (dendrite) shown in Fig. [T^, D = (dy) is the matrix corresponding 
to the nearest-neighbor coupling (cf. (2.26)). For branched dendrites (see Fig. [jj>), D may have a more 



complex structure. In either case, the structure of D reflects the geometry of the neuron. By combining this 



information, we obtain a model in the form of ( 1.2 1, (2.1 1, with 



11) - - Joi' 
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0, the coupling 
admits an alternative 



where 6\ = g^ 1 ^. The coupling is full rank. If we disregard calcium diffusion, i.e., set 5 
becomes partial with L = diag (1, 0). System of equations (4.1 ) and (4.2) with 6 

interpretation. One can view v^> and vy> as the membrane potential and calcium concentration of Cell i 
in a neuronal population. The coupling is due to the current through the gap-junctions between adjacent 
cells. In this case, D reflects network connectivity. If gap-junctional conductance depends on voltage or 
calcium concentration, the coupling is nonlinear as in (3.1 1. If the gap-junctions permit ions in one direction, 
coupling matrix D is not symmetric. 



In the remainder of this section, we present the results of numerical simulations of (4.1 ) and (4.2 1. We 
choose the variant of the model with partial coupling, i.e., with 6 = 0. When 5 > 0, the model has even 
better synchronization properties. The upper panel of Fig. [2] shows the phase plane and the time series of 
five uncoupled oscillators forced by small noise. The initial condition is chosen on the limit cycle of the 
deterministic system and is the same for each oscillator. Fig. [2^ shows phase trajectories of all five oscillators 
for approximately one cycle. As one can see from Fig. [2J5, after a few first cycles, under the influence of 
noise the oscillations gradually loose coherence. In contrast, the trajectories shown in Fig. |2fc and d remain 
tightly bundled. In these simulations, we used nearest-neighbor D2 and all-to-all coupling Di respectively. 

To illustrate Theorem[2] for simulations shown in Figure[2£, we use a random dissipative matrix. Specif- 
ically, we pick the entries of (— Q) 1//2 from a uniform distribution on [0, 1]: 



(-Q) 



1/2 



/ 





7577 





7060 





8235 





4387 





4898 \ 







7431 





0318 





6948 





3816 





4456 







3922 





2769 





3171 





7655 





6463 







6555 





0462 





9502 





7952 





7094 


V 





1712 





0971 





0344 





1869 





7547 / 



and let 



D fl 



QAo 



V 



-1 


0251 


2 


2043 


-1.6032 


0.5044 


-0.0804 


-0 


1264 





2772 


-0.3006 


0.2060 


-0.0562 


-1 


1549 


2 


5819 


-1.9613 


0.5210 


0.0133 


-0 


8807 


1 


9231 


-1.0823 


0.0333 


0.0066 


-0 


9049 


1 


8778 


-1.0060 


0.3772 


-0.3441 



(4.8) 



The trajectories in Figure|2^ are not as close to each other as in the two previous plots. To see how the value 
of k computed for different network topologies correlates with the degree of coherence in our simulations, 
we compute k(Di) = 0.8, k(D2) = 2, and k(Ds) = 23.1675. In accord with (2.34 1, the numerics show 
that coherence is better for smaller values of k. 

In conclusion, we relate the class of dissipative matrices to that of matrices that have been previously 
known to promote synchrony. The coupling matrices analyzed in j9j[T5l are subject to the constraint that the 
off diagonal elements are nonnegative. Note that many of the off diagonal elements of our randomly picked 
dissipative matrix D3 are negative. A straightforward albeit tedious calculation shows that if one chooses 
a dissipative matrix at random in the way we did in this example, the probability that at least one (or for 
that matter any fixed) off diagonal element is negative, is positive. This shows that the class of dissipative 
matrices is substantially bigger than those satisfying sufficient conditions for synchronization in ||9l [T31 . 



5 Discussion 



In this Letter, we have identified dissipative operators, a class of linear coupling operators that enforce syn- 
chrony in networks of oscillators provided that the interactions between oscillators are sufficiently strong. 
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Figure 2: Numerical simulations of compartmental model (4.1 )-(4.2). (a) Phase trajectories of five uncou- 
pled oscillators are plotted for approximately one cycle, (b) Timeseries of five uncoupled oscillators. All 
oscillators were given identical initial condition lying on the limit cycle of the deterministic system. Even 
small noise leads to desynchronization of five oscillators already after first several cycles. In contrast, simu- 
lations of the system coupled using nearest-neighbor (c), all-to-all (d), and random dissipative (e) coupling 
matrices show good coherence. The coherence is better for smaller values of k. 
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Our results apply to a broad class of networks including those with asymmetric, time-dependent, and nonlin- 
ear separable coupling schemes; as well as networks of local systems with nonperiodic attractors. Further- 



more, we have derived an analytic estimate ( 2.37 ) for the coherence of the network dynamics in the presence 
of noise. Robustness to noise is one of the main indicators gauging physical feasibility of the dynamical 
regimes generated by mathematical models. In this respect, ( 2.37 > gives important practical information 
about the factors contributing to the robustness of synchronous oscillations to noise. On the other hand, 
stability of the relevant dynamical states is among the key parameters determining the asymptotic value of 



the variance of the trajectories of a randomly perturbed dynamical system. For large systems like ( 1.2 1, an- 
alytical estimates of the quantities characterizing stability (e.g., Lyapunov exponents) are rare. By studying 



the variability of the synchronous regime in a randomly perturbed problem ( 1 .2 1 and (2.1), one can infer the 
degree of stability of the synchronous solution of the underlying deterministic system. Specifically, smaller 



values of /c(D) imply better stability of the synchronous solutions of ( 1.2 )o and (2.1 1. Importantly, k(D) 
reveals the contribution of the network topology to the stability of the synchronous state. Therefore, using 



the randomly perturbed model ( 1.2 1, (2.1 1 and the main estimate (2.37 1 may be viewed as a probabilistic 
method for studying stability of the synchronous solutions in the deterministic system. 
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Appendix. Parameter values for (4.1) and (4.2) 



The equations for the local systems in the neural network (4.1 1 and (4.2 1 are adopted from a nondimensional 
model of a dopamine neuron ll22l . For biophysical background and details of nondimesionalization, we 
refer an interested reader to ll22l . For the purposes of the present Letter, the values of several parameters 
of the original model were modified to make the oscillations less stiff. The parameter values used in the 
simulations shown in Figure [2] are summarized in the following table. 

Table 





1 


E 2 


-0.9 


£3 


-0.3 


9i 


0.8 


92 


2 


93 


1 


9 


0.3 


ai 


-0.35 


a-2 


1.4- 10" 2 


03 


1.8 


e 


0.1 


T 


5.0 




5.0 


a 


0.0001 
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